Computing the current density j in the coil

\(\def\curl{\operatorname{curl}}\def\Curl{\operatorname{Curl}}\def\div{\operatorname{div}}\) Before we solve the non-linear magnetostatic problem, we first need to find the current \(j\) flowing inside the coil \(\Omega_c\).

The first property is that \(j\) is solenoidal, i.e. \(\div(j)=0\). Further, we presume the existence of a potential \(\phi\) with \(j=-\sigma\nabla\phi\), where \(\sigma\) denotes the connectivity. As for the boundary conditions, we prescribe \(n\cdot j = -\sigma\partial_n\phi=0\) on the exterior boundary \(\Gamma_{ex}\), and \(\phi=1\) and \(\phi=0\) on the inflow \(\Gamma_{in}\) and outflow boundary \(\Gamma_{out}\), respectively. Altogether, we have

\[\begin{split}\begin{align} j + \sigma\nabla\phi &= 0 \qquad\text{on }\Omega_c \\ \div j &=0 \qquad\text{on }\Omega_c \\ n\cdot j = \sigma\partial_n\phi &= 0 \qquad\text{on }\Gamma_{ex} \\ \phi &= 0 \qquad\text{on }\Gamma_{out} \\ \phi &= 1 \qquad\text{on }\Gamma_{in} \end{align}\end{split}\]

However, in our case, the coil is a loop with no inflow or outflow boundary. This is where the face “coil_cut_1” we defined in the geometry comes into play. For this purpose, we have to introduce “fictitious” points and introduce a clone of the face “coil_cut_1” in order to be able to prescribe the necessary boundary conditions.

Ok so how do we even solve this? Weak formulation leads to

\[\begin{align} -\Delta(\sigma\phi)=0 \end{align}\]

The boundary condition \(\partial_n\phi=0\) on \(\Gamma_{ex}\) is natural, while the Dirichlet conditions are essential and have to included in the space. Further, since we solve the problem only on the coil, the other DOFs have to be eliminated from the final system.

We begin by loading the geometry and the mesh generated in the previous document

%%capture
%run TEAM_13_geometry.ipynb

Create the MESH object from the mesh generated by netgen

import sys
sys.path.insert(0,'../../../') # adds parent directory
import pde
MESH = pde.mesh3.netgen(geoOCCmesh)
MESH

np:4132, nt:23369, nf:3059, ne:720, nf_all:46966, ne_all:27728

The piece of code below duplicates the face, as described above, and generates a new MESH object.

import numpy as np

face_index = pde.tools.getIndices(MESH.regions_2d,'coil_cut_1')
faces = MESH.f[MESH.BoundaryFaces_Region == face_index,:3]
new_faces = faces.copy()

points_to_duplicate = np.unique(faces.ravel())
new_points = np.arange(MESH.np, MESH.np+points_to_duplicate.size)

actual_points = MESH.p[points_to_duplicate,:]

t_new = MESH.t[:,:4].copy()
f_new = MESH.f[:,:3].copy()
p_new = MESH.p.copy()

for i,pnt in enumerate(points_to_duplicate):

    # append point to list
    p_new = np.vstack([p_new,p_new[pnt,:]])

    # finding tets coordinates containing the ith point to duplicate
    tets_containing_points = np.argwhere(t_new[:,:4]==pnt)[:,0]

    for _,j in enumerate(tets_containing_points):
        #check if tet is left
        if MESH.mp_tet[j,0]<0:
            t_new[j,t_new[j,:]==pnt] = MESH.np + i

    # finding faces containing the points
    faces_containing_points = np.argwhere(f_new[:,:3]==pnt)[:,0]
    for _,j in enumerate(faces_containing_points):
        #check if face is left
        if 1/3*(p_new[f_new[j,0],0] + p_new[f_new[j,1],0] + p_new[f_new[j,2],0])<0:
            f_new[j,f_new[j,:]==pnt] = MESH.np + i


    # print(faces_containing_points)

t_new = np.c_[t_new,MESH.t[:,4]]

for i,j in enumerate(faces.ravel()):
    new_faces.ravel()[i] = new_points[points_to_duplicate==j][0]

new_faces = np.c_[new_faces,np.tile(MESH.f[:,3].max()+1,(new_faces.shape[0],1))]
f_new = (np.r_[np.c_[f_new,MESH.f[:,3]],new_faces]).astype(int)

regions_2d_new = MESH.regions_2d.copy()
regions_2d_new.append('new')

identifications = (np.c_[points_to_duplicate,new_points]).astype(int)
# stop
MESH = pde.mesh3(p_new,MESH.e,f_new,t_new,MESH.regions_3d,regions_2d_new,MESH.regions_1d,identifications = identifications)
MESH

np:4149, nt:23369, nf:3076, ne:720, nf_all:47028, ne_all:27806

As we can see, additional points and faces have been created to accommodate the additional interface.

Next, we proceed to solve the problem in \(H^1\). In weak form, we have:

Find \(\phi\in V^* = \{u\in H^1\;:\; u|_{\Gamma_{out}}=0 \text{ and } u|_{\Gamma_{in}}=1\}\) such that \((\sigma\nabla\phi,\nabla v) = 0\), for all \(v\in H^1\)

After homogenization, which involves splitting the solution \(\phi = \phi_* + \phi_0\), where \(\phi_*\in V^*\), we can instead solve:

Find \(\phi_0\in H^1\setminus\Gamma_{in,out}\) such that \((\sigma\nabla\phi_0,\nabla v) = -(\sigma\nabla\phi_*,\nabla v)\), for all \(v\in H^1\setminus\Gamma_{in,out}\)

After using conforming \(P_1\) continuous finite elements, we solve the discretized system in the following way: we split the discrete solution \(\phi = R^T_{in}\phi_{in} + R^T_{out}\phi_{out} + R^T_{int}\phi_{int} = R^T_{in}\phi_{in} + R^T_{int}\phi_{int}\), where \(R^T_*\) matrices which assign the correct indices to the boundaries and the interior degrees of freedom. Assume the stiffness matrix is denoted by \(K\) and the right-hand side is \(r\). Then, purely in \(H^1\), the system is of the form \(K\phi=f\). Plugging in the splitting for \(\phi\) leads to the system \(KR_{int}^T\phi_{int} = -KR_{in}^T\phi_{in}\). By multiplying from the left with \(R_{int}\), we obtain the square system

\[\begin{align} R_{int}KR_{int}^T\phi_{int} = -R_{int}KR_{in}^T\phi_{in} \end{align}\]

order = 1
D = pde.int.assemble3(MESH, order = order)
DB = pde.int.assembleB3(MESH, order = order)
N1,N2,N3 = pde.int.assembleN3(MESH, order = order)
unit_coil = pde.int.evaluate3(MESH, order = order, coeff = lambda x,y,z : 1+0*x, regions = 'coil')

###########################################################################

phi_H1 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'M', order = order)
dphix_H1, dphiy_H1, dphiz_H1 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'K', order = order)
phiB_H1 = pde.h1.assembleB3(MESH, space = 'P1', matrix = 'M', shape = phi_H1.shape, order = order)

R0, RS0 = pde.h1.assembleR3(MESH, space = 'P1', faces = 'new,coil_cut_1')
R1, RS1 = pde.h1.assembleR3(MESH, space = 'P1', faces = 'coil_cut_1')

M = phi_H1 @ D @ unit_coil @ phi_H1.T

K = dphix_H1 @ D @ unit_coil @ dphix_H1.T +\
    dphiy_H1 @ D @ unit_coil @ dphiy_H1.T +\
    dphiz_H1 @ D @ unit_coil @ dphiz_H1.T

r = -RS0 @ K @ R1.T @ (1+np.zeros(R1.shape[0]))
K = RS0 @ K @ RS0.T


RZ = pde.tools.removeZeros(K)
K = RZ @ K @ RZ.T
r = RZ @ r

sigma = 1#58.7e6
from sksparse.cholmod import cholesky as chol
p = chol(sigma*K).solve_A(r)
potential_H1 = RS0.T @ RZ.T @ p + R1.T @ (1+np.zeros(R1.shape[0]))

jx_L2 = -(dphix_H1.T@potential_H1)*unit_coil.diagonal()
jy_L2 = -(dphiy_H1.T@potential_H1)*unit_coil.diagonal()
jz_L2 = -(dphiz_H1.T@potential_H1)*unit_coil.diagonal()

dphix_H1_P0, dphiy_H1_P0, dphiz_H1_P0 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'K', order = 0)
unit_coil_P0 = pde.int.evaluate3(MESH, order = 0, coeff = lambda x,y,z : 1+0*x, regions = 'coil')
jx_L2_P0 = -(dphix_H1_P0.T@potential_H1)*unit_coil_P0.diagonal()
jy_L2_P0 = -(dphiy_H1_P0.T@potential_H1)*unit_coil_P0.diagonal()
jz_L2_P0 = -(dphiz_H1_P0.T@potential_H1)*unit_coil_P0.diagonal()

Instead of solving in \(H^1\), we can instead solve the mixed problem directly by employing \(H(\div)\) conforming finite elements.

order = 1
phix_Hdiv, phiy_Hdiv, phiz_Hdiv = pde.hdiv.assemble3(MESH, space = 'RT0', matrix = 'M', order = order)
divphi_Hdiv = pde.hdiv.assemble3(MESH, space = 'RT0', matrix = 'K', order = order)

phi_L2 = pde.l2.assemble3(MESH, space = 'P0', matrix = 'M', order = order)

D = pde.int.assemble3(MESH, order = order)

M_Hdiv_coil_full = phix_Hdiv @ D @ unit_coil @ phix_Hdiv.T +\
                   phiy_Hdiv @ D @ unit_coil @ phiy_Hdiv.T +\
                   phiz_Hdiv @ D @ unit_coil @ phiz_Hdiv.T

C_Hdiv_L2 = divphi_Hdiv @ D @ unit_coil @ phi_L2.T
R1, RS1 = pde.hdiv.assembleR3(MESH, space = 'RT0', faces = 'coil_face')

M_Hdiv_coil_full = RS1 @ M_Hdiv_coil_full @RS1.T
C_Hdiv_L2 = RS1 @ C_Hdiv_L2

import scipy.sparse as sp
AA = sp.bmat([[M_Hdiv_coil_full, C_Hdiv_L2],
              [C_Hdiv_L2.T, None]])

RZdiv = pde.tools.removeZeros(AA)
AA = RZdiv @ AA @ RZdiv.T

phiB_Hdiv = pde.hdiv.assembleB3(MESH, space = 'RT0', matrix = 'phi', shape = phix_Hdiv.shape, order = order)
unit_coil_B = pde.int.evaluateB3(MESH, order = order, coeff = lambda x,y,z : 1+0*x, faces = 'new').diagonal()
DB = pde.int.assembleB3(MESH, order = order)

rhs = unit_coil_B @ DB @ phiB_Hdiv.T
rhs = np.r_[RS1@rhs,np.zeros(MESH.nt)]
rhs = RZdiv @ rhs

xx = sp.linalg.spsolve(AA,rhs)

potential_L2 = (RZdiv.T@xx)[-MESH.nt:]
j_hdiv = RS1.T@(RZdiv.T@xx)[:-MESH.nt]

jx_hdiv = (phix_Hdiv.T@j_hdiv)*unit_coil.diagonal()
jy_hdiv = (phiy_Hdiv.T@j_hdiv)*unit_coil.diagonal()
jz_hdiv = (phiz_Hdiv.T@j_hdiv)*unit_coil.diagonal()

##############################################################################
unit_coil_P0 = pde.int.evaluate3(MESH, order = 0, coeff = lambda x,y,z : 1+0*x, regions = 'coil')

phix_Hdiv_P0, phiy_Hdiv_P0, phiz_Hdiv_P0 = pde.hdiv.assemble3(MESH, space = 'RT0', matrix = 'M', order = 0)
jx_hdiv_P0 = (phix_Hdiv_P0.T@j_hdiv)*unit_coil_P0.diagonal()
jy_hdiv_P0 = (phiy_Hdiv_P0.T@j_hdiv)*unit_coil_P0.diagonal()
jz_hdiv_P0 = (phiz_Hdiv_P0.T@j_hdiv)*unit_coil_P0.diagonal()

##############################################################################

grid = pde.tools.vtklib.createVTK(MESH)
pde.tools.vtklib.add_H1_Scalar(grid, potential_H1, 'potential_H1')
pde.tools.vtklib.add_L2_Scalar(grid, potential_L2, 'potential_L2')
pde.tools.vtklib.add_L2_Vector(grid,jx_L2_P0,jy_L2_P0,jz_L2_P0,'J_L2')
pde.tools.vtklib.add_L2_Vector(grid,jx_hdiv_P0,jy_hdiv_P0,jz_hdiv_P0,'J_HDIV')
pde.tools.vtklib.writeVTK(grid, 'current_density.vtu')

import pyvista as pv
mesh = pv.read('current_density.vtu')
mesh

Header	Data Arrays
UnstructuredGrid Information N Cells 23369 N Points 4149 X Bounds -2.000e+02, 2.000e+02 Y Bounds -2.000e+02, 2.000e+02 Z Bounds -1.000e+02, 1.000e+02 N Arrays 5	Name Field Type N Comp Min Max potential_H1 Points float32 1 0.000e+00 1.000e+00 Scalars_ Cells float64 1 0.000e+00 5.000e+00 potential_L2 Cells float32 1 0.000e+00 9.960e-01 J_L2 Cells float32 3 -2.110e-03 2.123e-03 J_HDIV Cells float32 3 -1.902e-03 1.931e-03

mesh.set_active_scalars("Scalars_")
threshed = mesh.threshold([0,1])

p = pv.Plotter()
threshed.set_active_scalars("potential_H1")
p.add_mesh(threshed, style='surface', opacity = 0.4, label=None)

mesh.set_active_vectors("J_L2")
arrows = mesh.glyph(scale="J_L2", orient=True, tolerance=0.03, factor=9500.0)
p.add_mesh(arrows, color="black")

p.camera_position = [(0, 0, 600),(0, 0, 0),(0, 0, 0)]
p.show(jupyter_backend='html')

mesh.set_active_scalars("Scalars_")
threshed = mesh.threshold([0,1])

p = pv.Plotter()
threshed.set_active_scalars("potential_L2")
p.add_mesh(threshed, style='surface', opacity = 0.4, label=None)

mesh.set_active_vectors("J_HDIV")
arrows = mesh.glyph(scale="J_HDIV", orient=True, tolerance=0.03, factor=9500.0)
p.add_mesh(arrows, color="black")

p.camera_position = [(0, 0, 600),(0, 0, 0),(0, 0, 0)]
p.show(jupyter_backend='html')